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The propagation of lower hybrid waves in the presence of ponderomotive density fluctu- 
ations is considered. The problem is treated in two dimensions and, in order to be able 
to correctly impose the boundary conditions, the waves are allowed to evolve in time. 
The fields are described by iv T — J <iC + u CC + \ v \ 2 v = where v is proportional to the 
electric field, r to time, and C and £ measure distances across and along the lower hybrid 
ray. The behavior of the waves is investigated numerically. If the amplitude of the waves 
is large enough, the spectrum of the waves broadens and their parallel wavelength 
becomes shorter. The assumptions made in the formulation preclude the application of 
these results to the lower hybrid heating experiment on Alcator-A. Nevertheless, there 
are indications that the physics embodied in this problem are responsible for some of 
the results of that experiment. 
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I. INTRODUCTION. 

The injection of rf power near the lower hybrid frequency is an attractive method for 
the auxiliary heating of tokamak plasmas. 1 Because of the high powers required (several 
MW) and because lower hybrid waves principally propagate along well-defined resonance 
cones, there has been considerable interest in nonlinear effects on the propagation of lower 
hybrid waves. This problem was first addressed by Morales and Lee 5 who studied the two- 
dimensional electrostatic propagation of one of the two lower hybrid rays in a homogeneous 
plasma. Although this is perhaps the simplest model that can be considered, a correct 
treatment of this problem has yet to be made. It is this deficiency that this paper attempts 
to remedy. 

Briefly the difficulty of this problem arises as follows: If it is assumed that the rf fields in 
the plasma have reached a steady state, i.e., that the potential is given by <f>(x, z) exp(— iuj^t) 
(x and z are coordinates perpendicular and parallel to the ambient magnetic field Bo), then 
the electric field obeys the complex modified Korteweg-deVries equation. 5 ~ 7 This equation is 
mathematically well-posed when solved as an initial value problem in one of the coordinates 
x, that is when (f>(x = 0, z) is given. Unfortunately, this does not correspond to physically 
realizable boundary conditions since waves can propagate in both the +x and — x directions 
in a single ray. 7 ' 8 When the correct boundary conditions are imposed, there is numerical 
evidence that solutions of the complex modified Korteweg-deVries equation need not exist 
and that this equation is therefore ill-posed. 7 This is confirmed by our finding solutions 
inconsistent with the assumption of a steady state (see Sec. IV). The problem arises because 
the direction of power flow which determines how to impose the boundary conditions is 
defined only with reference to a problem in which a temporal evolution of the wave packet is 
allowed. In assuming a steady state for the electric field amplitude, the equation no longer 
has built into it the crucial ingredient which determines how the boundary conditions are 
imposed. This defect is corrected by including a slow time dependence of the potential [so 
that the potential is given by 4>(x, z, t) exp(— iuiot)]. This then leads to a nonlinear partial 
differential equation in two spatial dimensions and time, which we will study numerically in 
this paper. 

In formulating this problem, we shall ignore many effects which should possibly be 
included to obtain a complete understanding of the propagation of lower hybrid waves. 
This will enable us to study the effect of the nonlinearity in as simple a system as possible. 
Even so, the numerical solution of a partial differential equation in two dimensions and 
time is time consuming and, as we shall see, the behavior of the fields as described by this 
equation can be quite complicated. In the absence of any analytical methods for solving 
this equation, we must therefore be content with the solution for only a few cases. This 
will enable us to confirm the threshold for strong nonlinear effects given in Ref. 7 and to 
give a more complete description of the nature of these nonlinear effects. Because of the 
approximations made in formulating the problem, we shall not be able to apply these results 
to the propagation of lower hybrid waves near the edge of a tokamak plasma where nonlinear 
effects are most strong. However, there are indications that the physical processes that are 
considered in this paper do play a role in the propagation in that region. 

The plan of this paper is as follows: In Sec. II we derive the partial differential equation 
governing the temporal evolution of <f>. The basic properties of this equation will be discussed 
in Sec. III. In particular, we will show how the right boundary conditions automatically 
drop out from the equation. Since the equation is analytically intractable, we resort to a 
numerical integration, the results of which are given in Sec. IV. The threshold condition 
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for a strong nonlinear interaction is derived in Sec. V. The results are summarized and the 
consequences to lower hybrid heating of tokamaks are presented in Sec. VI. 



II. FORMULATION OF THE PROBLEM. 

In this section, we derive the partial differential equation governing the evolution of the 
electric field for one of the lower hybrid rays in two spatial dimensions and time. We assume 
that the fields are electrostatic, that the plasma is homogeneous, and that it is immersed 
in a uniform magnetic field. The derivation closely follows that of the complex modified 
Korteweg-deVries equation by Morales and Lee; 5 the additional ingredient we consider is 
the slow time dependence of 4>. 

The potential of the lower hybrid wave is taken to have the form 

Re[0(x, z, t) exp(— ioj^t)] 

where the t dependence of <j) is taken to be much slower than luq. Since we are only doing the 
two-dimensional problem, we have d/dy = 0. In the electrostatic limit, 4> obeys Poisson's 
law 

V-K(V,0/0t, |V0| 2 ) • = 0. (1) 

Here we regard the dielectric tensor K as an operator through its arguments V and d/dt. 
The dependence on |V0| accounts for the nonlincarity. If the temporal evolution of 4> 
is sufficiently slow (compared with ion acoustic time scales), then nonlinearities due to 
ponderomotive density changes can be written in this way. Parametric instabilities are 
excluded from our consideration. 

If the dependence of <f> on x, z, and t is weak and if the 4> itself is small so that the 
nonlinearity is weak, we may expand K to obtain 

K(V,3M,|V0| 2 ) = K + i^^:VV 



+ e ^7!77^ + l^l 2 + 0(£ 2 ). (2) 



<9K d <9K 

d{d/dtjdi + e d\V^f 

All the evaluations of K on the right hand side of Eq. (2) are at the point V = d/dt = 
|V<^| 2 = 0. The term involving <9K/<9V is zero because of the symmetries of a stationary 
plasma. We have introduced a formal expansion parameter e to aid in the ordering of terms. 
The fact that the last three terms in Eq. (2) are taken to be of order e results a maximal 
ordering. If, in fact, one of the terms is much smaller than the others, a subsidiary ordering 
can be introduced to eliminate that term. 

Since it is more usual to write K as a normal function, rather than as an operator, 
we rewrite K as K(k, u, n) where k is the Fourier-transform variable conjugate to space, 
u> is that conjugate to time, and n is the plasma density. [A space-time dependence of 
exp(ik • r — iwt), where r = (x, z), is assumed.] Then the derivatives in Eq. (2) become 

<9 2 K d 2 K dK _ .dK dK dn dK 

~ ~<9k<9k' d{d/dt) " * a? <9|V0| 2 ~ <9|V</>| 2 ^' 

In this case, all the evaluations on the right hand sides are performed at k = 0, w = u>o, 
and n — no, where no is the unperturbed density of the plasma (i.e., in the absence of any 
electric fields). 
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We now write down Eq. (1) to various orders in e. To order e°, we find 

V- K(0,wo,no)- V0 = O (3) 

where 

K(0,w ,no) = . 

Assuming that K±K\\ is negative, Eq. (3) is hyperbolic and its solution is 

4> = 4> + (x- gz) + (f>-{x + gz) 

where g = {-K^/ K^) 1 ' 2 . 

In general, both the right- and left-going rays will be present. However, if the source 
of the lower hybrid rays is localized, the two rays will be separated far from the source and 
in that region we may then treat each ray in isolation. Similarly, only a single ray will be 
present if the source is adjusted so that only that ray is excited. We may therefore assume 
that, to zeroth order, only the right-going ray is present, i.e., 4>- = 0. To order e , wc 
write (j) = <j>(et', ex', z') where t' = t, x' = x, z' = cz — sx, and s = (1 + g 2 )^ 1 ^ 2 , c = gs. 
We should think of x' and z' as measuring distances along and across the lower hybrid ray. 
Substituting this form of </> into Eq. (1), we obtain to order e 1 

, d 2 E dE d 3 E d „ l2 N 

where E = d<j)/dz' is the electric field measured in the direction normal to the lower hybrid 
ray, A, B, C, and D are real coefficients given by 

8K 

A = k — k, B =2sK ± , 

ouj 

C=-ikk:5i:kk, D 1 9n 



2 ' Okdk' ' ' " n9|V0| 2 ' 

and k is the unit vector (— s, c). 

These expressions may be evaluated using the expression for K for a stationary Maxwell- 
ian plasma. 9 Usually the frequency of the incident rf power satisfies fif <§; 0Jq <C fl 2 where 
Clj is the cyclotron frequency of species j (j = i or e for ions or electrons). In that case, the 
components of the zeroth-order dielectric tensor may be written 

K ± = l + ^-^4, K n = 1 - ^ e Wpi 



where uj p j is the plasma frequency for species j with density no. The coefficients in Eq. (4) 
become 




S C —hr—— + 3c 



An {T & + Tiy 
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where vfj = Tj/mj, Tj is the temperature of species j, and eo is the dielectric constant 
of free space. In writing the expression for D, we have used the formula for the density 
depression due to the ponderomotive force. 5 ' 10, 11 Note that C is always positive. The fact 
that C and D have the same signs has an important effect on the propagation of the lower 
hybrid waves. 

We now rescale the variables in Eq. (4): 

Z = x'/(BVC), r = t'/A, ( = z'/Vc, v = VDE. (5) 
The equation for v is then 

f C 2 

iv T — / d( + + \v\ u = (6a) 

J — oo 

where subscripted Greek letters denote differentiation. When integrating Eq. (4) to give Eq. 
(6a), we assume that v and all its derivatives are zero at £ = — oo (i.e., far from the ray). 
The accessibility condition 12 imposes an auxiliary condition on v namely 



f 

J — C 



vd( = 0; (6b) 



in other words, the electric potential is equal at £ = ±oo. Equation (6) describes the 
evolution of the electric field of a single lower hybrid ray in two dimensions and time under 
the influence of thermal dispersion and nonlinear ponderomotive effects. Equation (6a) is 
the same as Eq. (50) of Rcf. 7 if the replacements v — > v* , a — > r, r — > £, and £ — > — £ are 
made to the latter equation. 

In order to give some appreciation of the scaling of the variables in Eq. (5), we give 
more explicit forms for them in the limit uj^ C w§ < Wp e <C O 2 an d ^ (these may 
be the conditions applicable as the lower hybrid wave propagates through the low density 
part of a tokamak plasma). Equation (5) then becomes 

p _ x ujpt gz-x _( e \ 1/2 

where X De = v te /uj pe is the Debye length and g = w / w pe- 



III. BASIC PROPERTIES OF THE EQUATION. 

Equation (6) is invariant under the scaling transformation 

£^A~ 3 £, t^\- 2 t, C^A^C, v^Xv. (8) 

This invariance allows us, for example, to normalize the scale length of the ray in the £ 
direction to unity. 

If we set d/dr = in Eq. (6), we obtain, as expected, the complex modified Korteweg- 
deVries equation as the equation obeyed by steady-state fields. If <9/(9£ = 0, we obtain 
the nonlinear Schrodinger equation, which is soluble by the inverse scattering method. 13 
Unfortunately, not much use can be made of this fact because, as we shall see, the imposition 
of the correct boundary conditions involves solving Eq. (6) in a finite domain in £. 

For each of the four conservation laws obtained for the complex modified Korteweg- 
deVries equation, there is a corresponding conservation law for Eq. (6). These are given in 
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nable I. The second and third of these laws are statements of the conservation of momentum 
and energy respectively. We will return to these shortly. 

In order to determine the correct boundary conditions for Eq. (6), we Fourier transform 
in £. Defining the Fourier transform of v by 



V(£,k,t) = u(£,C,r)exp(-wcC)d< 

J — oo 

we obtain 

V T + cV e + iQV + iN{V) = 0, (9a) 
V(k = 0) = 0, (9b) 

where c = 1/k, il = k 2 , and N(V) is the Fourier transform of — \v\ 2 v. This equation is 
hyperbolic in (r, £) space. However, the direction of the characteristic velocity c depends on 
k. This means that we must specify boundary conditions at each end of a strip in £. More 
precisely, the full initial and boundary conditions required for solving Eq. (9) in the domain 

£o < £ < £l, — OO < £ < OO, T <T<OO 

are 

F(£o <£<£i,k,t ), V(£ ,k>0,t>t ), y(a,«<0,r>r o ). (10) 

Taking the second and fourth conservation laws in Table I, integrating them over £ and 
£, and transforming them to k space, we obtain the laws of conservation of momentum and 
energy as they apply to the domain of the problem 



d /-£i r°° r 00 £i 

— / / Mdndi+ / Tdn 
dr J£ J-oo J-oo 



= 0, (11) 

£ = £o 
£i 

= 0, (12) 

£ = £o 



where M = k\V\ 2 is the spectral momentum density, T — cM is the spectral force density, 
£ = \V\ 2 is the spectral energy density, and V = c£ is the spectral power density. The 
characteristic velocity c is the group velocity of the lower hybrid waves. The boundary 
conditions given in Eq. (10) stipulate that we should specify the waves entering the domain 
of the problem, but not those that leave this domain. 

Given that we are treating a nonlinear problem, it may be a little surprising the bound- 
ary conditions are applied in the same way as for the corresponding linear problem. The 
mathematical reason for this is that the nonlinear term does not involve any r or £ deriva- 
tives. This arises because, in deriving Eq. (6), we assumed that the nonlinear term was very 
small. The expressions for the energy density, power density, etc. are therefore correctly 
given by linear theory. Although the nonlinearity is small, it can have a finite effect because 
the ordering also assumes that the length over which the interaction takes place is large. 

Notice that the nature of the equation has forced us to take boundary conditions which 
are in accord with physical notions of power flow. This has happened because we have the 
time variation in Eq. (6). If we consider the steady-state equation, the complex modified 
Korteweg-deVries equation, which is obtained by setting djdr = in Eq. (6) 

v ( - v ccc - (\v\ 2 v) c = 0, (13) 
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then we would clearly like to impose the same boundary conditions as for Eq. (6) namely 
Eq. (10). But the concept of power flow cannot be defined from Eq. (13) and so there is no 
guarantee that the boundary conditions given in Eq. (f 0) are sensible. Indeed, the evidence 
of Rcf. 7 suggests that Eq. (13) is ill-posed with these boundary conditions. This can be so 
because, even if Eq. (6) is given steady-state boundary conditions, i.e., 

d d 
< £ < f i , k, to -> -co), g^V^o, k > 0, r) = , « < 0, r) = 0, 

it is not necessarily the case that the solution to Eq. (6) approaches a steady state. 

As mentioned in Sec. II, the ordering we used in deriving Eq. (6) was a maximal 
ordering. We can eliminate some of the terms in Eq. (6) using a subsidiary ordering. E.g., 
the linear limit is obtained by letting v — > Sv and taking the limit 5 — > 0. Similarly, the 
non-dispersive limit is obtained by letting £ — > £/<5, £ — > £/<5 and again taking the limit 
6 — > 0. It is instructive to examine this case in more detail because there are indications 
that if £i — £o is sufficiently large then the dispersive term must be included. Consider the 
steady-state equation (13) without the dispersive term v^. Multiplying this equation by 
v* and adding the complex conjugate (these are the same operations that lead to the second 
conservation law in Table I), we obtain 

U£ + uuq = (14) 

where u — 3|u| 2 . Since the boundary conditions involve argu, Eq. (14) should be supple- 
mented by another equation. However, this is not so if we impose initial conditions v(£ = 
£o, 0- As is well known, the solution develops a shock at £ = £ + 1/ max[c?w(£ = £ , 
where du/dC, is infinite. Before this point is reached the ordering assumed for the inde- 
pendent variables, i.e., that d/d( = 0(d), breaks down. Inclusion of the dispersive term 
prevents the shock from forming. Indeed, Eq. (13) when treated as an initial value problem 
in £ is well-posed in an infinite domain in £. Although Eq. (6) with the correct boundary 
conditions is different from Eq. (13) with initial conditions, it is likely that Eq. (6) is also 
ill-posed if the dispersive term is ignored and if the width of the domain in £ is large enough. 



IV. NUMERICAL RESULTS. 

We investigate the solution to Eq. (6) by numerically solving Eq. (9) with boundary 
and initial conditions 

V(0<£ <&,t = 0) = 0, (15a) 
V(£ = 0, k > 0, t > 0) = v a TTu(n - Ko)(k — Ko) cxp[— j(k — Ko) 2 ], (15 b) 

V(£ = £i,k<0,t>0) = 0, (15 c) 

where u is the unit step function, vo and Ko are real constants and kq > 0. The method of 
solution is given in Appendix A. Without loss of generality, we have taken t = £o = 0. In 
£ space Eq. (15b) becomes 

«(£ = °)Uo = «0[1 + C z (0\ exp(moC) 

where Z is the plasma dispersion function. With this form of the boundary conditions, we 
now have three parameters to vary: v 0} Ko, and A£ = £i — £o = £i- The geometry of the 
waveguide array which might give this input spectrum will be discussed in Sec. V. 
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We begin by considering the linear non-dispersive limit; i.e., we remove the third and 
fourth terms in Eq. (6a). Figure 1 shows the solution at various times. Early in the evolution 
of the waves, only the lowest k components have propagated throughout the plasma because 
these have the greatest group velocities. Later on, the other k components have had time 
to reach the far end of the system, and the steady state is attained. As t — > oo, we 
have v = v(Q. In the original (x,z) coordinate system, this corresponds to the familiar 
propagation along resonance cones. 

In order to further diagnose the results, we identify three sets of waves, the incident 
waves (£ = and k > 0), the reflected waves (£ = and k < 0), and the transmitted 
waves (£ = £i and k > 0). With our choice of boundary conditions the incident waves are 
constant for r > and there are no waves incident at the boundary £ = £i. We define 
reflection and transmission coefficients, R and T, as the ratios of the instantaneous power 
reflected and the instantaneous power transmitted to the incident power. These powers are 
defined by integrating V over the appropriate segment of the boundary. Note that, because 
energy can be stored in the plasma, we need not have R + T = 1 , although this is true 
in a time-averaged sense. Furthermore, it is possible that R or T can momentarily exceed 
unity. In order to determine where in k space the power is concentrated, we define (k) for 
each wave component as the ratio of the force to the power for that component. (Since 

2 2 

V = \V\ /k and T = \V\ , this ratio has the dimensions of k.) 

For the case shown in Fig. 1, R is obviously zero and the average k of the reflected 
waves, (/c) r , is undefined. Figure 2 shows T and the average n of the transmitted waves, 
(«)t, for this case. The transmission coefficient T rises in steps as each k mode reaches the 
far boundary. [Because of the use of periodic boundary conditions in the numerical solution 
of Eq. (6), there are only a discrete set of k modes.] Similarly, {n) t rises with time since 
modes with higher values of n take longer to traverse the domain. As t — > oo, T approaches 
unity, signifying total transmission, and (n) t approaches (ir/2) 1 / 2 w 1.25 which is the same 
as the average value of k for the incident waves, («;),. 

We now turn to the nonlinear dispersive problem. Because we expect that the effect of 
the nonlinearity will increase with both vo and A£, we begin with a case where these quan- 
tities are chosen small enough that the nonlinearity only weakly modifies the propagation. 
Figure 3 shows the solution for such a case with vo — 3, ko — 0, and A£ = 0.1. It is found 
that v attains a steady state by r = 2 approximately and only this state is shown in Fig. 
3. The evolution to this steady state is similar to that for the linear equation. In Fig. 4 
are displayed the reflection and transmission coefficients. Note that the reflection coefficient 
settles down to a small value (somewhat less than 2%). 

As the amplitude of the boundary value vo is increased, the solution exhibits quite 
different behavior from the linear solution. In Fig. 5, we have chosen the same parameters 
as in Fig. 3, except that vq has been increased to 4. Now v appears never to reach a 
steady state, but instead oscillates in some aperiodic fashion. This can be seen in Fig. 
6, where R, T, («) r , and (n)t are plotted. From Fig. 5, we can see what is happening 
during these oscillations. Various k components in the forward wave nonlinearly interact to 
produce a reflected wave [Fig. 5(b)]. This interaction causes a severe depletion of the low n 
components of the forward wave and transfer of this energy into higher k components of the 
forward wave and into the reflected wave [Fig. 5(c)]. After the interaction has nearly gone 
to completion, the nonlinearly excited components of the field transit out of the system 
and the fields relax to a state in which there is little variation in |V| with £ [Fig. 5(d)]. 
The nonlinear interaction then begins again and the cycle approximately repeats itself [Fig. 
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5(e)]. In £ space, this interaction is manifested by a narrowing and peaking of the electric 
field amplitude v. A typical field pattern is shown in Fig. 5(f). 

Returning to Fig. 6, we see that the reflection is appreciable, oscillating around about 
20%. Correspondingly, T oscillates around 80%. Because the energy stored in the field 
comes out in bursts, T is occasionally greater than unity. Finally, we note that {n)t exceeds 
{n)i s=a 1.25 by approximately a factor of 2. 

Because the long-time solution of the fields does not reach a steady state, it is useful to 
describe the solution in terms of its temporal spectrum. To do this, we take a time record of 
the transmitted and reflected waves for r a < r < r b . We choose r a to be large enough for the 
transients to have disappeared from the solution and rt — r a to be large enough to include 
several oscillations of the solution. Effects arising from the finite record length — r a are 
partially eliminated by multiplying by a cosine window i{l — cos[27r(r — T a )/(T& — r a )]}. 
The resulting function is then transformed in time using a discrete Fourier transform to give 
the output spectrum V(£' , k, a) where a is the slow frequency variable (conjugate with r), 
£' = £i for k > (the transmitted wave), and £' = for n < (the reflected wave). [We 
assume a space-time dependence of exp(in( — iar).] 

In Fig. 7, we plot the output spectrum V(£', k, cr) for the case shown in Figs. 5 and 6. 
The spectrum is computed between r a — 1 and t& = 5. The large peaks at cr = are the 
steady-state components of the transmitted and reflected waves. The broad spectrum in a 
is symptomatic of the aperiodic nature of the waves. There is an interesting correlation in 
the spectrum: the positive frequency (a) components of both the transmitted and reflected 
waves tend to have larger values of \k\ than the negative frequency components. 

We have seen that the solution becomes more turbulent and that the reflection increases 
as «o is increased. The other parameters that may be varied are A£ and Ko- We now 
investigate the influence of these parameters. 

We begin by varying A£. Figure 8 shows the output spectrum for vq = 4, ko = and 
A£ = 0.05. We see that there is almost no energy in the components of the wave with 
cr 7^ 0; i.e., the fields nearly reach a steady state. The reflection coefficient is smaller than in 
the previous case (where A£ = 0.1) being approximately 4%. Similarly {n) t is only slightly 
greater than (1.4 as compared to 1.25). On the other hand, if we increase the value 
of A£ to 0.2 (twice its value in Fig. 7), a broad turbulent spectrum is recovered (Fig. 9). 
Although R and {n) t have approximately the same values as for A£ = 0.1, there is now 
more energy in the non-steady components of the output spectrum. 

Figures 10 and 11 show the effect on the output spectrum of increasing k (keeping 
v Q = 4 and A^ = 0.1). In Fig. 10, where k = 1, we see that spectrum consists of a few 
narrow bands indicating that the solution is nearly periodic. Furthermore, the reflected 
wave consists almost entirely of negative frequency components. The reflection coefficient 
for this case is about 12%, or half of what it was with k = 0. The increase of («;) t over 
is moderate (about 3 as compared to 2.42). Increasing k to 2 (Fig. 11), we see much the 
same picture except that the frequency of the oscillations is higher. The reflection coefficient 
of is further reduced to 7% and there is now little difference between (n) t and («;},. 
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V. THRESHOLD FOR THE NONLINEARITY. 

In the case shown in Figs. 5-7, the nonlincarity causes three phenomena of importance 
to lower hybrid heating experiments: (1) the reflection coefficient is appreciable; (2) the 
solution reaches a turbulent state; (3) the mean wavenumber of the waves transmitted into 
the plasma is increased. Comparing the results presented here with those of Sec. VIII of 
Ref. 7 in which the solution of the steady-state problem was attempted, we see that the 
conditions for the occurrence of these phenomena closely agree with the conditions under 
which the steady-state problem had no solution and reflection was large. 

Before quoting the threshold results from Ref. 7, we use the scaling invariance Eq. (8) 
to rewrite the boundary condition Eq. (15b) as 

V(£ = 0, k > 0,t > 0) = v nu(n- ko)(k-ko) AC 2 exp[-±AC 2 (K - k ) 2 ]. 

In real space, we have 

»(£ = °)Uo = Ml + (C/AC)£(C/AC)] cxp(z«oC). 

This roughly corresponds to excitation of lower hybrid waves by a waveguide array of width 
A(. The amplitude of the electric field in the plasma is v and the phasing of the waveguides 
is such as to give a wavenumber of kq. Thus, if the waveguides are phased 0, tt, 0, tt, . . . , 
we have kq = ttM/A( where M is the number of waveguides. [The average wavenumber 
of the spectrum of a single lower hybrid ray is (n)i which, for M large, is approximately 
ko + 2(2/7r) 1//2 /A£. The first term is attributable to the phasing of the waveguides and the 
second term arises because of the finite width of the waveguide array.] 

For these boundary conditions, the threshold conditions given in Ref. 7 are 

v > 2ko, (16a) 
AC > V2w,7 3 . (16b) 

The latter conditions is only correct for M small (< 4). For larger values of M, it should 
be replaced by 

A£>ACV 2 , (16c) 

which specifics that A£ should exceed the distance for a shock to form in the dispersion- 
less equation (14). Equation (16) gives the conditions under which the three phenomena 
described at the beginning of this section occur. 

If we now undo the normalizations using Eq. (7), we find that Eq. (16) becomes 

E z0 Sz/T e > 4V3tt, (17a) 

Ax > ±V6(T e /E x0 )/p, (17b) 

Ax>gAz/l3, (17c) 

where E x0 and E z0 are the electric field amplitudes perpendicular and parallel to the mag- 
netic field (so that E z q = gE x o), Az is the width of the waveguide array, 5z — Az/M is 
the width of a single waveguide (assuming a 0, n, 0, n, . . . phasing), Ax is the width of the 
nonlinear region of the plasma in the x direction, and (3 is j€oE x0 /n T e , the ratio of the 
electric field energy to the plasma kinetic energy. 
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VI. DISCUSSION. 

We have examined the nonlinear evolution of a lower hybrid wave in two dimensions 
and time. Under the conditions given in Eq. (17), the nonlinearity can cause appreciable 
reflection, turbulent variation in the fields, and an increase in the wavenumbers of the 
transmitted waves. In Ref. 7, it was found that these conditions could be satisfied in a 
small laboratory plasma or near the edge of a tokamak plasma. In such circumstances, the 
assumption of a steady state and the analyses in Refs. 5-7 based on this assumption are 
wrong. 

The most important application of lower hybrid waves is for heating a tokamak plasma 
to ignition temperatures. It is therefore important to understand the processes that may 
modify the lower hybrid waves before they have penetrated to the center of the plasma. 
Unfortunately, it is not possible to obtain quantitative estimates of these effects with the 
theory as outlined in the preceding sections. The reason is that, as we have pointed out, the 
nonlinearity is only important in a tokamak near the edge of the plasma and in that region 
many other physical processes are likely to be involved. Examples of effects which should 
probably be included to give a full understanding of lower hybrid wave propagation near the 
edge arc: nonlinear coupling to the second (left-going) ray; electromagnetic effects; density 
and temperature gradients; saturation of the nonlinearity; ion inertia in the low frequency 
equations; coupling to low-frequency drift waves. 

Although Eq. (6) should be regarded as a very simplified model equation describing 
the propagation of lower hybrid waves near the edge of a tokamak plasma, some of the 
phenomena predicted by this equation have been observed in the lower hybrid heating 
experiment on Alcator-A. 14 ' 15 The C02-laser scattering diagnostic on that experiment 
indicated a broadening of the frequency spectrum of the waves. This is consistent with 
the turbulent spectra seen in the solutions of Eq. (6). Furthermore, the spectrum becomes 
asymmetric as the wave propagates into the plasma, the peak of the spectrum being shifted 
down from the frequency of the injected waves. This may be a result of the observation in 
Sec. IV that the components of the transmitted and reflected waves which are down-shifted 
in frequency have a smaller (in absolute value) wavenumber than those which are up-shifted. 
If those waves with higher wavenumbers are preferentially damped as the wave travels into 
the plasma (e.g., by electron Landau damping), wc would expect to see a net downwards 
shift in the spectrum of waves. 

The other somewhat puzzling result of the Alcator-A experiment was the apparent 
independence of the heating results to the phasing of the waveguides. In addition, in order 
to explain the density dependence of the heating results, it was postulated that the n z 
(parallel index of refraction k z c/uj) spectrum of the waves was shifted to being peaked at 
around n z — 5 irrespective of the phasing of the waveguides. (This is to be compared with 
values of n z predicted on the basis of linear theory of less than 1.5 for the waveguides in phase 
and of 3 for the waveguides out of phase.) Again, this is qualitatively in agreement with 
the theoretical results of Sec. IV. There, the average wavenumber of the waves transmitted 
from the nonlinear region into the center of the plasma was roughly twice that of the wave 
injected into the nonlinear region when kq = (corresponding to the wave guides being in 
phase). As the phase difference between neighboring waveguide becomes finite (i.e., as n 
is increased), the amount by which the wavenumber spectrum is shifted is reduced. 

It is not clear whether the nonlinear reflection predicted by the theory given in this 
paper would be observed as an increase in the reflected power measured in the waveguides. 
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It may be that this power is reflected again on the cutoff at u> = w pe very close to the plasma 
edge. This would convert the power into the other lower hybrid ray. 

Thus, it appears that the physics included in Eq. (6) may be responsible for some of the 
results of the Alcator-A experiment. In order to be able to say definitively whether or not 
the experimental observations are a result of the processes included in Eq. (6), two advances 
are required. Firstly, the theory has to be refined to a point where quantitative comparisons 
are possible. Secondly, additional experimental observations are required. (At present, 
Alcator-A is the only tokamak experiment for which the detailed observations provided by 
the C02-laser scattering diagnostic are available. In addition, it would be helpful to have 
accurate measurements of the plasma conditions near the edge of the plasma.) 

We should point out that other theories have beeen advanced to explain the Alcator-A 
results. In particular, Bonoli and Ott 16 have suggested a linear theory. This is supported 
by the observed linear dependence between the density fluctuations and the applied power 
(over a fairly large range) which suggests either that the phenomena are linear or else that 
the nonlinear processes saturate at a low amplitude. 

Another interesting theoretical result, given by Morales, 17 concerns the coupling of rf 
energy into lower hybrid waves at the plasma edge. A density gradient was included in his 
model and a temporal evolution was allowed. As in the theory presented here, it was found 
that a steady state need not be reached; rather, the rf energy entered the plasma in bursts 
similar to the behavior seen in Fig. 5. However, only a single n z component was included 
so that the nonlinear coupling of different n z components was disallowed. 

In summary, we have presented a theory of the nonlinear propagation of lower hybrid 
waves. At sufficiently large powers, the fields become turbulent and the wavenumber spec- 
trum is shifted upwards. The theory should accurately describe the propagation in small 
laboratory devices. While we may expect qualitatively similar results near the edge of a 
tokamak plasma, other physical effects need to be included to obtain an accurate description 
of the lower hybrid fields in this region. 
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APPENDIX A: NUMERICAL PROCEDURE FOR SOLVING EQ. (6). 

The numerical solution of Eq. (6) is carried out in Fourier space so that Eq. (6) becomes 
Eq. (9). Periodic boundary conditions are used in the £ direction, v(£, ( + L,t) = v(£, (, r); 
therefore, the Fourier spectrum is discrete the spacing between modes being Sk = 2tt/L. 
We use n points to describe v over a single period L. Transformations between the k and 
£ spaces are achieved using a discrete Fourier transform. Thus, V is approximated by n 
Fourier modes. The £ coordinate is approximated by a grid whose spacing is 5£. This 
equation is solved by splitting each time step St into two equal pieces and by approximating 
Eq. (9a) in the interval < t < St by 

fcVf, for < t < Ut, , . 

V = — 2 x < ~~ (Al) 

\iQV + iN(V), for \5t <t < St. v ' 

(For simplicity, we describe the solution only for the first time step. The extension beyond 
this is obvious.) 

In each half time step this equation is a partial differential equation in only two inde- 
pendent variables. During the first half time step, k is merely a parameter and we have a 
simple linear wave equation to solve. We approximate V = V(t = \St) as given by Eq. 
(Al) by shifting V(t = 0) over C = cSt/S^ grid positions. This step is exact if C is an 
integer. Normally, this is not the case and in that event we interpolate between neighboring 
grid points. We use linear interpolation on the quantities \ V\ 2 and \ V\ 2 argV. This ensures 
that in this step energy and momentum are conserved. The reason for interpolating in 
|F| 2 argF rather than a,rgV is that in the former case the ambiguity of the argument of 
V = is irrelevant. The grid positions within C of the boundary are set to the boundary 
value. 

During the second half time step we have to solve the nonlinear Schrodinger equation 
at each position £. The method is similar to that used in Ref. 7 to solve the steady-state 
equation (13). The dispersive term is treated exactly and the nonlinear term is treated with 
a second-order Runge-Kutta scheme. Thus, V(t — St) is approximated by 

V" = BV' + DN(V'), 

V(5t) = BV + \D[N(V") + N{V% 

where B = exp(— iViSt) and D = —(1 — B)/fl. The nonlinear term ^V(T^) is calculated 
by transforming V into v using the discrete Fourier transform, computing — \v\ v, and 
transforming back into k space. In order to avoid problems of aliasing in the computation 
of N(V), the highest n/2 Fourier modes are artificially set to zero. 

The accuracy of the numerical integration is checked using the momentum- and energy- 
conservation laws. Specifically, the time integrals of Eqs. (11) and (12) are numerically 
computed. These are divided by the total momentum input into the plasma (at £ = £o) an d 
by the total energy injected into the plasma (both at £ = £ , « > and at £ = £i, k < 0) to 
provide two measures of the accuracy of the numerical integration, AM and AS. 

The results presented in Sec. IV were (with the exception of Figs. 1 and 2) computed 
with n = 2 8 , L = 20, and J£ = 10~ 3 . The time step was taken to be St = 1CT 3 in Figs. 
5-7, and 9, St = 2 x 1CT 3 in Figs. 8, 10, and 11, and St = 5 x 10~ 3 in Figs. 3 and 4. The 
error parameters AM and A£ were less than about 5 x 10~ 4 at t = 5 in all these cases. 

The case shown in Figs. 5-7 was also computed with n — 2 9 , L — 20, and St = <5£ = 
5 x 10~ 4 (i.e., the grid spacing was halved in r, £, and £). The solution agreed well with 
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the solution obtained using the coarser grid until about r = 1.5. Thereafter, the solutions 
diverged from each other. This is to be expected in a system which exhibits turbulent 
solutions because the solution is typically very sensitive to the initial conditions. Numerical 
errors, which have an effect similar to changing the initial conditions slightly, can therefore 
lead to large changes in the solution. However, although the detailed solution is different 
after r = 1.5 in these two cases, the general character of the solution is the same. Thus, 
the output spectrum for the case with the finer grid spacing shows the same features as the 
spectrum given in Fig. 7. 
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Tables 

TABLE I. Conservation laws for Eq. (6). The form of the conservation laws is dTj/dr + 
dXj/d£ + dZj/d( — 0. Here q is given by — v. 



j Tj Xj Zj 

1 ivq —v + I -17 1 2 v 

2 -iv*V£ \v\ 2 3 \v<;\ 2 - |w 2 | (( , - | \v\ 4 + ivv* 

3 iv*v/: \v c \ 2 - \ \v\ 4 - iv*v T -2Rc(v*v c ) + \q ( \ 2 

4 \v\ 2 —iq*v 21m(v* v^) + iq*q^ 




FIG. 1. The solution to Eq. (6) with the last two terms omitted. Here v = 1, k — 0, and 
A£ = 1. The solution is shown in both the k and ( spaces for [(a) and (b)] r = 0, [(c) and 
(d)] r = 0.5, and [(e) and (f)] r = 5. 




FIG. 2. The reflection coefficient R (a) and mean value of n transmitted {n) t (b) as func- 
tions of time for the case shown in Fig. 1. 




FIG. 3. The solution to Eq. (6) for w = 3, n n = 0, and A£ = 0.1. The solution is shown 
for t = 5 in (a) n space and (b) Q space. 



19- 




«0 



o.oto - 



0.005 - 



* 1 ' 1 1 — 1 i i < i i 

0.5 t.O 1.5 2.0 2.5 3.0 3.5 4.0 4.5 5.0 

r 



1.0 
0.8 
0.6 
0.4 
0.2 



(b) 



1 " 1 L 



■1 1 I 1 1 1 l 

0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0 4.5 5 

T 



FIG. 4. The reflection (a) and transmission (b) coefficients, R and T, as functions of time 
for the case shown in Fig. 3. 
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FIG. 8. The output spectrum for v = 4, k = 0, and A£ = 0.05. 




FIG. 10. The output spectrum for vq — 4, kq — 1, and A£ = 0.1. 



FIG. 11. The output spectrum for v — 4, k — 2, and A£ = 0.1. 



